Conservation  properties  for  the  Galerkin  and  stabilised  forms 
of  the  advection-diffusion  and  incompressible  Navier-Stokes 

equations 


Thomas  J.R.  Hughes1  *  Garth  N.  Wells2 


A  common  criticism  of  continuous  Galerkin  finite  element  methods  is  their 
perceived  lack  of  conservation.  This  may  in  fact  be  true  for  incompressible 
flows  when  advective,  rather  than  conservative,  weak  forms  are  employed. 
However,  advective  forms  are  often  preferred  on  grounds  of  accuracy  despite 
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ily  remedied,  and  conservative  procedures  for  advective  forms  can  be  devel¬ 
oped  from  multiscale  concepts.  As  a  result,  conservative  stabilised  finite  ele¬ 
ment  procedures  are  presented  for  the  advection-diffusion  and  incompressible 
Navier-Stokes  equations. 
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]  Introduction 


Conservation  is  a  highly  sort  after  property  in  numerical  simulations  of  transport  phe¬ 
nomena.  Global  conservation  is  intuitively  appealing  and  often  used  for  validating  a 
computational  model.  A  common  criticism  of  continuous  Galerkin  methods  is  that  they 
are  not  globally  or  locally  conservative.  This  misunderstanding  arises  from  the  problem 
that  it  is  not  possible  in  the  general  case  to  set  the  weight  function  equal  to  unity  in  a  con¬ 
tinuous  Galerkin  method  as  the  weight  function  must  vanish  where  Dirichlet  boundary 
conditions  are  applied.  Recently,  this  issue  has  been  clarified  by  Hughes  et  al.  [1],  where 
it  is  was  shown  that  continuous  Galerkin  methods  for  advection-diffusion  equations  are 
both  globally  and  locally  conservative. 

For  advection-diffusion  problems,  when  the  weak  form  is  written  in  advective,  rather 
than  conservative,  form  and  the  flow  field  is  divergence-free  (at  least  in  an  appropriate 
weak  sense),  global  conservation  of  the  advected  quantity  is  preserved.  However,  this 
condition  is  not  often  satisfied  by  stable  finite  element  methods  based  on  the  Galerkin 
formulation  and  when  the  velocity  field  is  computed  from  a  pressure-stabilised  finite  ele¬ 
ment  formulation.  This  leads  to  a  lack  of  conservation  of  the  transported  quantity,  despite 
the  method  being  convergent. 

It  is  shown  here  how  global  conservation  can  be  preserved  when  using  common  sta¬ 
bilised  finite  element  methods  through  a  small  residual-based  modification  of  existing 
formulations.  As  a  starting  point,  global  conservation  for  stabilised  methods  applied  to 
the  advection-diffusion  equation  is  proven  for  the  case  where  the  advective  velocity  is 
known  to  be  solenoidal.  The  examination  is  then  extended  to  the  case  where  the  velocity 
comes  from  the  solution  of  a  stabilised  incompressible  flow  problem  and  the  weak  form 
is  in  the  advective,  rather  than  conservative,  form.  It  is  shown  that  satisfaction  of  the 
inf-sup  condition  will,  in  general,  preclude  conservation  for  classical  Galerkin  methods. 
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and  stabilisation  of  the  continuity  equation  in  common  stabilised  finite  element  formu¬ 
lations  results  in  global  conservation  not  being  guaranteed.  Considering  the  governing 
equations  in  a  multiscale  context  [2,  3]  points  the  way  to  a  formulation  which  is  globally 
conservative.  The  procedure  is  then  generalised  to  the  Navier-Stokes  equations.  It  is  again 
shown  that  the  addition  of  particular  residual-based  terms  to  the  stabilised  formulation 
leads  to  correct  conservation  of  momentum.  Implementation  of  the  procedures  requires 
only  minor  additions  to  existing  stabilised  finite  element  codes. 

2  Conservation  for  the  scalar  advection-diffusion  equation 
2. 1  Strong  and  weak  forms 

Before  examining  conservation  properties,  it  is  useful  to  define  carefully  the  problem  (for 
background,  see  Hughes  et  al.  [4]).  Consider  O,  an  open,  bounded  region  in  \Rd,  where  d 
is  the  spatial  dimension.  The  boundary  of  Q  is  denoted  by  T  =  30,  and  the  outward  unit 
normal  to  T  is  denoted  n  =  . . .  ,?n).  Given  a  solenoidal  advective  velocity  field  a 

(V  •  a  =  0),  consider  the  following  definitions: 


(1) 

(2) 

(3) 


Note  that  a+  is  equal  to  a„  at  an  outflow  boundary,  and  is  equal  to  zero  at  an  inflow 
boundary.  Conversely,  a, 7  is  equal  to  an  at  an  inflow  boundary,  and  is  equal  to  zero  at  an 
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Figure  1:  Definition  of  boundary  partitions, 
outflow  boundary  Let  { F  ,  T+  }  and  { f  /, }  be  partitions  of  Y,  defined  by: 


r  =  {x  g  r  |  an  ( x )  <  0} 
r+  =  r  -  r~ 


inflow  boundary 
outflow  boundary 


(4) 

(5) 


and 


Tf  =  Ygn  T± 

r±  =  r,nr±. 


(6) 

(7) 


The  various  partitions  are  illustrated  in  Figure  2.1. 

Denoting  the  diffusivity  as  k,  which  is  assumed  to  be  constant  and  positive  (k  >  0),  the 
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following  fluxes  are  defined: 


cra  ( <p )  =  —cup 

(8) 

<xd  (c p )  =  KS/(p 

(9) 

( 

7  =  aa  +  ad 

(10) 

o n  =  n  ■  cra 

(11) 

crd  =  n  ■  ad 

(12) 

o 

n  =  n  ■  a. 

(13) 

Advective-diffusive  transport  of  a  scalar  cp  is  assumed  to  be  governed  by: 

1 

< 

II 

in  Q 

(14) 

<P=S 

on  Tg 

(15) 

-an(p  +  oi  (<p)  =  h 

on  T/; 

(16) 

where  /  :  Q  — >  R,  g  :  Tg  — >  R  and  h  :  T;,  - 

—>  R  are  prescribed.  The  boundary  condition  on 

F/,  can  be  interpreted  as: 

/ 

h  =  < 

h  on  T;; 

l 

(17) 

(  h+  on  T+ 

where  h  is  the  total  flux  and  h+  is  the  diffusive  flux. 


an  (<p)  =  h 

total  flux  boundary  condition 

(18) 

oil  (?)  =  h+ 

diffusive  flux  boundary  condition. 

(19) 

For  the  variational  formulation  of  the  advection-diffusion  equation,  the  following  func- 
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tion  spaces  are  required: 


5  =  j<£  G  H1  (Q)  \<p  =  g  on  rg} 

(20) 

V  =  |  we  H 1  (Q)  re  =  0  on  . 

(21) 

It  is  useful  to  define  another  function  space  W,  such  that 

H1  (fi)  =  V  ©  W 

(22) 

The  standard  variational  problem  for  advection-diffusion  consists  of:  find  (p  G  S  such  that 

B  (re,  (p)  =  L  (re)  V  re  G  V 

(23) 

where 

B  (re, (p)  =  (Vre,rr(^))n+  (zv,a+cp)Yh 

(24) 

and 

L  (re)  =  (w,f)n  +  (w,h)Th. 

(25) 

Consistency  of  the  variational  form  is  easily  proven  by  integrating  equation  (23)  by 
parts. 


0  =  B  (re,  u)  —L  (re) 

=  -  (re,  V  ■  cr(cp)  +  /)n  +  (re,  -a~  +  a*  (<p)  -  h) 

\  /  lh 


(26) 


which,  under  suitable  smoothness  hypotheses,  implies  the  original  strong  form  of  the 
governing  boundary-value  problem.  Stability  can  be  proven  by  showing  that  B  (iv,  (p)  is 
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coercive. 


B(w,w)  =  (Vzv,  —  azv  +  KVw)n  +  (w,a£zv)T 

1  1 

=  5  O,  (V  •  a)  zv) n  -  -  (i w,anw)Th  +  jc||  Vw||n  +  (ro,a+a;) 


=  -  (zv,  (V  •  a)  w)q  +  k||  Vai||n  +  -r 


\an\1/2zv 


(27) 


If  the  advective  field  is  divergence-free,  equation  (27)  is  equivalent  to: 


B  (w,w)  =  k\\ Vw||q  + 


\an\1/2iv 


2 

r„ 


(28) 


which  proves  that  the  stability  is  assured.  The  boundary  conditions  considered  here  were 
described  in  Hughes  et  al.  [4],  and  they  result  in  a  well-posed  variational  problem  if  the 
flow  field  is  divergence-free. 

Proving  conservation  for  the  weak  form  of  the  advective-diffusion  equation  in  both 
conservative  and  advective  formats  is  straightforward  for  the  case  Tg  =  0.  Considering 
equation  (23)  and  setting  zv  =  1, 


0  =  B  (l,(p)  —  L  (1) 

=  f  a£<p  dT  —  I  f  dd.  —  [  h  dT, 
Jr  in  Jr 


which  can  also  be  expressed  as: 


(29) 


0  =  -  (y  ( -ancp  +  h+ )  dT  +  J  h~  dT  +  J  f  dnj  (30) 

which  proves  that  <p  is  conserved  (recalling  that  h  '  is  the  diffusive  flux,  hence  (7/  1  —  an(p )  |r+ 
is  the  total  flux,  and  7/  is  also  the  total  flux).  Establishing  conservation  for  the  case  /  0 
is  more  complicated  as  it  is  not  possible  to  set  zv  =  1  in  the  entire  domain  since  zv  must  be 
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equal  to  zero  on  kg.  Proving  conservation  requires  the  definition  of  an  auxiliary  flux  h ,, 
on  the  boundary  Tg  [1],  Consider  the  problem:  find  cp  G  6"  and  hg  €  W  such  that 


B  (zv,cp)  =  L  (zv)  +  (zv,hg)T  V  zv  G  H1  (Cl) .  (31) 

Recalling  equation  (22),  this  problem  can  be  alternatively  expressed  as:  find  <p  G  S  and 
hg  €  W  such  that 


B  ( w ,  cp)  =  L  (zv)  V  zv  £  V 

(32) 

B(q,(p)  =  L(q)  +  (q,hg)Tg 

(33) 

Setting  zv  =  1  in  equation  (31),  which  is  possible  since  1  G  H1  (Q), 


0  =  B  (l,(p)  —  L  (1)  —  (l,/?g)r^ 

=  -{  [  ( -a„(p  +  h+ )  dY+  [  h~dT+  [  fdfl+  [  hg  dTj 
\Jt+  Jt-  J  n  -/r,  / 


(34) 


which  proves  that  cp  is  conserved. 

Conservation  has  been  shown  above  for  the  weak  form  of  the  advection-diffusion  equa¬ 
tion  in  conservation  form,  with  no  reliance  on  the  field  a  being  divergence-free.  From  the 
identity: 

V  •  (atp)  =  (V  •  a)(p  +  a  ■  V<p  (35) 


and  assuming  V  •  a  =  0,  the  weak  advection-diffusion  equation  in  advective  format  is 
expressed  as: 


Ba  (TV,<p)  =  ( zv,  a  ■  V0)n  +  (Vw,?cV0)n  -  (w,an<p)r  -  ( ivran(p)Yh  (36) 
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where  the  subscript  'a'  has  been  added  to  denote  the  advective  format.  Considering  the 
advective  form  and  setting  zv  =  1  for  the  case  Tg  /  0, 

0  =  Ba  (1  ,<p)  -  L  (1)  -  (zv,hg)  ry 

=  —  (1/  (V  •  a)  cp )n  +  —  (l,fl,;</>)r^  —  (l ,an  (P)yIj  ~  (l//)n  —  ~~  (l'^s)rg 

=  -  (1,  (V  •  «)  </>)n  +  -  (l//)n  -  (U)r*  - 

=  —  (  /  (V  •  a)  (p  dCl  +  /  (—an<p  +  fz+)  dT+  h~dT+  hg  dT  +  If  dCl  J  , 
yin  ./r+  Jrh  Jts  Jn  J 

(37) 

hence  cp  is  conserved  if  V  •  a  =  0. 

2.2  Conservation  in  Galerkin  and  stabilised  Galerkin  formulations 

Conservation  for  the  case  in  which  Dirichlet  boundary  conditions  are  applied  (T<s>  /  0)  in 
a  Galerkin  formulation  can  be  proven  in  a  similar  way  to  the  infinite-dimensional  case  in 
the  previous  section.  To  prove  conservation,  again  an  auxiliary  flux  must  be  defined  on 
Tg  [1],  Consider  now  Sh  C  S  and  Vh  C  V  to  be  finite  element  spaces.  The  space  Sh  con¬ 
tains  basis  functions,  each  associated  with  a  node,  which  satisfy  the  Dirichlet  boundary 
conditions  (technically  speaking,  Sh  is  a  linear  manifold,  not  a  linear  space).  The  space 
Vh  is  the  span  of  the  basis  functions  associated  with  each  node  in  Q  excluding  those  lying 
on  I  , .  A  complementary  space  Wh  is  defined  which  is  the  span  of  of  the  basis  functions 
associated  with  all  nodes  lying  on  Tg.  It  is  assumed  that  the  support  of  the  traces  of  func¬ 
tions  in  Wh  is  contained  in  Tg.  A  further  space  is  defined  as: 

Xh  =  Vh  ©  >V,!.  (38) 

Hence,  X1'  is  the  span  of  the  shape  functions  associated  with  all  nodes.  Note  that  Sh  C  Xh. 
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For  generality,  conservation  will  be  shown  for  the  case  of  stabilised  methods,  which  are 
a  superset  of  the  Galerkin  method.  The  domain  Q  is  subdivided  into  'elements'  Cle,  and 
the  domain  of  element  interiors  O  is  defined  as: 

n  =  [J  ne  (39) 

e 

which  does  not  include  the  element  boundaries. 

The  standard  Galerkin  problem  consists  of:  find  <ph  €  Sh  such  that 

B  (wh,(ph^j  =  L  ( wh )  V  wh  €  Vh.  (40) 

A  class  of  stabilised  Galerkin  methods  can  be  expressed  as: 

B  (wh,  (ph)  +  (Cawh,  t ar)  =  L  (wh)  V  wh  G  Vh  (41) 

in  which  r  is  the  residual  of  the  advection-diffusion  equation, 

r  =  a  ■  V(ph  -  kA cph  -  f,  (42) 

Tfl  is  a  stabilisation  parameter  and  Ca  is  an  operator  which  defines  the  stabilisation  method. 
For  the  streamlined-upwind  Petrov-Galerkin  method  (SUPG)  [5],  Cn  is  given  by: 

Cawh  =  a  ■  Vwh.  (43) 

For  the  Galerkin/ least-squares  method  (GLS)  [4], 

Cnwh  =  a  ■  Vwh  -  kAzv’1  (44) 
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and  for  the  multiscale  method  (MS)  [2,  6] 


Cawh  =  a  ■  Vzvh  +  kAhA  (45) 

Expressions  for  Ta  may  be  found  in  References  [7-11]. 

Proving  conservation  for  the  case  Pg  =  0  for  the  presented  stabilised  Galerkin  methods 
is  trivial  as  the  stabilisation  term  involves  the  gradient  of  zvh,  which  vanishes  when  setting 
zvh  =  1  everywhere  in  Q.  For  the  non-trivial  case  of  F,s,  /  0,  consider  the  following 
Galerkin  problem:  find  <ph  G  Sh  and  hg  G  Wh  such  that 

B  (wh, <ph^j  +  (Cazvh,  rar^j ^  =  L  (zvh"j  +  (wh, hfj ^  V  wh  G  Xh  (46) 

where  hg  is  the  flux  across  F?.  This  problem  can  be  split,  yielding:  find  (ph  G  Sh  and 
hg  G  Wh  such  that 

B  (wh,  <ph^j  +  ( Cazvh ,  % -ar)  —  L  (V)  =0  V  zvh  G  V;'  (47) 

B  (qh,4>h)  +  rflr)  ^  -  L  (V)  =  (qh,hf)^  V  qh  G  W,!.  (48) 

Equation  (47)  is  the  standard  stabilised  advection-diffusion  problem,  which  can  be  solved 
without  knowing  hg.  Once  <ph  has  been  computed  from  equation  (47),  it  can  be  inserted 
into  equation  (48)  to  calculate  the  auxiliary  flux,  which  is  merely  a  post-processing  proce¬ 
dure  (see  Hughes  et  al.  [1]  and  the  references  therein). 

Conservation  for  the  case  /  0  can  be  proven  by  setting  zvh  =  1  in  equation  (46)  (note 
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that  zvh  =1  €  Xh,  irrespective  of  the  boundary  conditions). 


o  =  B  (l,  <ph)  +  {Cal,  t flr)  A  -  L  (1)  -  (l,  hfj 

=  I  a+(phdT-  f  h  dT  —  [  f  dCl  —  [  hh  dT 
J  Th  Jrh  Jo  Jrg 

=  +  h+)  dT  +  f  h~dT  +  Jfdn  +  J  h'^dT 


(49) 


which  proves  that  Galerkin  and  the  considered  stabilised  Galerkin  finite  element  meth¬ 
ods  are  globally  conservative.  While  the  present  arguments  are  of  a  global  nature,  the 
same  procedure  can  be  extended  to  local  subdomains  consisting  of  individual  elements 
or  unions  of  elements,  thus  establishing  that  continuous  Galerkin  and  stabilised  Galerkin 
methods  are  also  locally  conservative  [1]. 

Thus  far,  it  has  been  assumed  that  a  is  solenoidal.  Consequently,  the  conservative  and 
advective  forms  of  the  weak  advection-diffusion  equations  are  interchangeable.  Often, 
the  advective  form  is  preferred  in  practice.  However,  a  problem  arises  in  the  advective 
form  when  a  is  computed  numerically  from  another  equation,  such  as  the  incompressible 
Navier-Stokes  equations,  and  fails  to  satisfy  the  divergence-free  condition.  In  this  case. 


(vzvb,  a(ph^j  n  ~  (jv11'  (V  •  a)(ph^J  +  ( 'zvh,an(ph ^  (50) 


in  which  the  term  involving  V  •  a  is  not  identically  zero.  For  conservation  to  hold,  the 
term  involving  V  •  a  in  equation  (50)  must  vanish  when  zvh  =  1.  That  is, 

[  (V  •  a)  (ph  dCl  =  0.  (51) 

Jo 

In  a  Galerkin  finite  element  method  for  the  incompressible  Navier-Stokes  equations,  this 
condition  will  be  satisfied  when  the  space  of  pressure  interpolations  Vh  contains  Xh,  be- 
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cause,  in  this  case,  the  discrete  approximation  of  incompressibility  is  given  by: 

f  qh  (V  •  a)  dd  =  0  V  qh  £  Vh  (52) 

Jn 

which  implies  that  equation  (51)  holds  as  Vh  Xh.  Unfortunately,  Vh  ^  Pi'1'  is  often 
the  case  due  to  restrictions  imposed  by  the  inf-sup  condition  (see  Brezzi  and  Fortin  [12]). 
Furthermore,  it  is  a  common  practice  to  use  some  form  of  pressure  stabilisation  in  which 
case  equation  (52)  definitely  will  not  hold.  In  order  to  construct  conservative  methods 
for  these  cases,  an  approach  based  on  multiscale  considerations  proves  fruitful,  and  is 
pursued  herein. 


2.3  Conservation  when  using  a  computed  flow  field  from  a  stabilised  Galerkin 
method 

To  construct  a  conservative  advection-diffusion  scheme  for  cases  when  the  advective  flow 
field  comes  from  a  numerical  simulation,  the  starting  point  is  a  multiscale  decomposition 
of  the  advective  flow  field.  Following  the  multiscale  concept  [2,  3],  the  velocity  field  is 
decomposed  additively, 

a  =  a  +  a'  (53) 

where  a  is  the  'coarse'  scale  component  of  the  velocity  field  and  a '  is  the  'fine'  scale  com¬ 
ponent  of  the  velocity  field.  It  is  assumed  that  a  =  a  on  T.  Returning  to  equation  (50)  and 
inserting  the  multiscale  decomposition. 


(w,  a  ■  V(p)n  =  (zv,  ( a  +  a')  ■  Vip)^ 

=  (w,V-  (a<p))n-  {zv,(V  ■  a)(p)n+  (w,a’ ■ 


(54) 
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Integrating  the  first  term  by  parts  on  the  right-hand  side  of  equation  (54),  and  considering 
that  a  =  a  on  T, 

(zv,  a  ■  V<^>)n  =  -  (Vzv,  d(p)n  +  ( w,an<p)T  —  (zv,  (V  •  a)  (p)n  +  (zv,  a'  •  V</>)n  .  (55) 

Inserting  this  expression  into  the  advective  form  of  the  advection-diffusion  equation  (36) 
and  including  the  boundary  conservation  term, 

0  =Ba  ( zv,<p )  -  L  (zv)  -  (■ zv,hg)T 

=  -  ( Viv,acp)a  +  ( w,a„(p)T  -  (zv,  (V  •  a)  < p )n  +  (zv,  a'  ■  V^)n  (56) 

+  (Vzv,K'V(p)Ci  -  ( zo,an<p)Tg  -  (zv,an(p)Th  -  ( zv,f)n  -  ( w,h)Th  -  (zv,hg) ^  . 

Setting  zv  =  1, 


0  =  (l,a+<p) Th  -  (1,  (V  •  a)  <p)a  +  {1,0!  ■  V*)n  -  (l,/)n  -  (1,/Or,  -  (57) 

which  can  be  equivalently  expressed  as: 


0  =  - 


0(— an(p  +  h+)  dT  +  [  fdCl+  I  h  dT  +  j  hg  dT 

J  o  « 

+  /  (V  •  «)  (j)  dCl  —  I  a'-'VzpdO, ).  (58) 

./n  J  o  / 


From  the  above  equation,  it  is  clear  that  (58)  reduces  to  (34),  and  (p  is  conserved,  if  the  last 
two  terms  cancel.  That  is,  if: 


((V  •  a)  <p  —  \7<p  ■  a ')  dCl 


=  0. 


(59) 
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Requiring  that: 


/  (q  (V  •  a)  -  Vq  ■  a')  dfi  =  0  \/q  £  H1  (Q)  (60) 

ensures  that  equation  (59)  holds,  and  conservation  is  guaranteed.  From  multiscale  con¬ 
siderations  and  in  light  of  the  connection  to  stabilised  finite  element  methods  [2,  3],  an 
approximation  of  the  fine  scale  velocity  on  element  interiors  is  given  by  [2]: 

a'  «  —rpr  (61) 

where  rp  is  a  stabilisation  parameter  (pressure  stabilisation)  and  r  is  the  residual  from 
the  equation  governing  the  advective  flow  field  (e.g.,  the  incompressible  Navier-Stokes 
equations).  Considering  now  that  a  is  represented  by  the  finite-dimensional  solution  a1', 
inserting  the  model  for  a '  into  equation  (60)  yields: 

(qh,  V  •  +  ( Vqh ,  rpr)  =  0  V  qh  G  Vh  (62) 

which  is  the  pressure-stabilised  continuity  equation  for  the  SUPG,  GLS  and  MS  methods. 
Hence,  a  conservative  form  for  advection-diffusion,  in  advection  format,  for  problems 
involving  an  advective  velocity  field  which  satisfies  equation  (62)  and  Vh  2  Sh,  is  given 
by: 

(wh,a  ■  +  (vwh,KV(ph^  -  (wh,an(ph^  -  (wh,a~(ph^ 

-  (V,  r pr  ■  Vcph^j n  +  (zvh, h'j r  +  (zvh, hfj r  Vivh  €  Xh.  (63) 

The  last  term  on  the  right-hand  side  of  equation  (63)  has  the  form  of  an  advection  term, 
hence,  in  addition  to  the  usual  advective  term,  it  too  needs  to  be  stabilised.  A  stabilised. 
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conservative  formulation  is  given  by: 


[zvh,a  ■  V(ph^j  ^  +  (vzvh,KV(ph^j  ^  ^ ivh,an(ph  j  ^ zvh,an(ph ^ 

+  [cnwh,rar^^  +  {zvh,a'  •  V<£,!)n  +  (C'awh,  r'  ( a '  •  V<£,!))n 

=  (V, /) ^  +  (whf h'j r  +  (wh, hfj r  Vivh  €  Xh  (64) 

where  a '  =  —  r^r,  r'  is  a  stabilisation  parameter  for  the  fine-scale  term,  and 


> 

V. 

1 

II 

> 

"<3 

II 

—  <3 

) 

(65) 

a  ■  Vzvh 

(SUPG) 

CgWh  =  <  a  •  \7ivh  _  KAlVh 

(GLS) 

(66) 

a  ■Vzvh  +  kAzvh 

(MS). 

The  first  term  in  the  box  is  the  key  to  conservation,  while  the  second  term  in  the  box 
stabilises  the  first  and  does  not  affect  conservation.  The  second  term  in  the  box  can  be 
equivalently  expressed  as: 

(£>,J,  r'  (V  •  V/)  )  ^  =  ( Cawh ,  n  (67) 

which  elucidates  its  symmetric,  positive-definite  nature.  Note  that  the  required  modifica¬ 
tions  to  an  existing  computer  code  for  stabilised  Galerkin  methods  to  ensure  conservation 
are  minor. 

Remark 

The  boxed  terms  are  consistent  modifications  of  the  method  due  to  the  fact  that  r  is  a 
residual  vanishing  for  the  exact  solution.  The  stabilisation  parameter,  t' ,  is  O  (h/\a'\)  = 
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O  (h/  ( Tp  | r | ) ) .  Consequently,  the  additional  term  vanishes  when  r  =  0.  A  more  precise 
expression  for  x'  is  presented  in  Taylor  et  al.  [13].  □ 

2.4  Time-dependent  case 

For  generality,  the  preceding  developments  are  extended  to  the  unsteady  case.  Consider 
the  following  space-time  domains: 


Q  =  Qx]0,T[  (68) 

Q  =  Qx]0,  T[  (69) 

P  =  rx]0,T[  (70) 

Pg  =  rgx]0,T[  (71) 

Ph  =  r;jx]0,T[  (72) 

P+  =  r+x]0,T[  (73) 

Ph~  =  r,;x]0,T[.  (74) 


Solution  of  the  unsteady  advection-diffusion  equation  involves:  find  <p  =  <p(x,t)  Vx  G 
O,  Vf  G  [0,  T]  such  that: 


^-V.<r(*)=/ 

in  Q 

(75) 

<p(x,  0)  =  (po 

in  Q 

(76) 

<P  =  g 

on  Pg 

(77) 

-a-(p  +  ai  (cp)  =  h 

on  Ph 

(78) 

where  /  :  Q  — >  1R,  g  :  Pg  — >  1R  and  /j  :  Pp,  — >  R  are  prescribed. 

For  the  discrete  problem,  the  space-time  domain  Q  is  divided  into  'time  slabs'.  Consider 
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the  following  space-time  slabs. 


Q"  =  Ox]f„,f„+1[  (79) 

Q"  =  Ox]f„,f„+i[  (80) 

p"  =  rx]f„,fn+1[  (81) 

Pl  =  Tgx]tn,tn+1[  (82) 

ph  =  rhx]tn,tn+i[  (83) 

P+"=r+x]fn/fn+i[  (84) 

Phn  =  T^x]tn,tn+1[  (85) 


where  0  =  to  <  t\  <  ■ . .  <  fw  =  T.  The  relevant  function  spaces  on  each  slab  are  similar 
to  those  for  the  steady  case,  with  infinite  dimensional  spaces  defined  by: 


Sn  =  {<f>sH1( Qn)  |  <^  =  y  on  P” |  (86) 

V”  =  jzc  G  H1  (Qn)  |  w  =  0  on  P”}  (87) 

H1  (Q'!)  =  V”  ©  Wn  (88) 

and  relevant  finite-dimensional  spaces  defined  by: 

(. Sn)h  C  Sn  (89) 

( Vn)h  C  V”  (90) 

(Sn)h  C  (. Xn)h  =  (Vn)h  ©  (Wn)h  .  (91) 
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Solution  of  the  unsteady  problem  requires:  find  cph  €  ( Sn)h  such  that. 


+  (a^,«'  •  V<^)Q„  +  (£>\r'  (a1  •  V^!))q„ 


"’*'/)  O.  +  K")p.  +  M)  ■*"  (In)) 


n 


Vw,!  g  (*”)*,  n  =  0, 1, . . . ,  N  -  1  (92) 


where  a'  =  —  rpr,  the  operator  C'a  is  defined  by  equation  (65),  and 


L-awh  =  < 


h 


dzv 

dt 

dwh 

dt 

dzvh 

l  dt 


+  a  ■  Vzvh  (SUPG) 

+  a  ■  Vzvh  —  KAivh  (GLS) 

+  a  ■  Vwh  +  kAiv h  (MS). 


(93) 


The  boxed  terms  are  the  same  as  for  the  steady-state  case  (see  equation  (64)).  Considering 


/  (V  •  a)  f1  dQ  +  /  x vr  ■  Vf  dQ  =  0, 

/  Q»  JQ” 


(94) 


which  holds  when  a  is  computed  from  a  pressure-stabilised  finite  element  method  and 
components  of  the  advective  velocity  field  a,-  come  from  the  same  space  as  zph,  and  setting 
zvh  =  1  in  equation  (92), 

[  <ph{t-+1)  dCl  =  f  <ph(t~)  dCt+  [  an<phdP  +  [  a~<phdP+  [  f  dQ 
Jn  Jo.  Jp’i  J  p"  J  Qn 

+  f  hdP+  [  hh  dP-  [  a-VzphdQ+  f  rpr  ■  V(ph  dQ.  (95) 
JP’’  JPn  "  J  Qn  J  Qn 
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Expanding  the  term  involving  a  and  considering  equation  (94), 


[  (ph{t~+1)  dd=  [  cf>h  (t~)  dn+  [  (-an<ph  +  h+)  dP 
J  n  J  n  J  p+"  v  ' 

+  [  h~dP+  [  fdQ+  f  ftJdP  (96) 

dP-"  JQ"  dp;; 

which  shows  that  cp1'  is  conserved  on  space-time  slabs. 

3  Conservation  for  the  incompressible  Navier-Stokes  equations 
3.1  Incompressible  Navier-Stokes  equations 

The  developments  in  this  section  for  the  incompressible  Navier-Stokes  equations  mirror 
the  developments  for  the  advection  diffusion  equation.  Before  proceeding  to  the  Navier- 
Stokes  equations,  it  is  useful  to  establish  some  definitions.  Denoting  the  velocity  field  u, 


un  =  n  ■  u 


(97) 


u 


+ 

n 


U 


n 


U-n  +  |  Un 

2 

Hjl  U  yi  | 

2 


(98) 

(99) 


where,  similar  to  before,  «+  is  equal  to  un  at  an  outflow  boundary  and  zero  elsewhere, 
and  u~  is  equal  to  u„  at  inflow  boundary  and  zero  elsewhere.  Let  {T“,r+}  and  {Tg,  T;,} 
be  partitions  of  T,  defined  by: 


T  =  {x  G  T  |  un  (x)  <  0}  inflow  boundary  (100) 

T+  =  F  —  T~  outflow  boundary  (101) 
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and,  as  before. 


Tf  =  r?  n  r±  (102) 

it  =  r/(  n  r±.  (103) 

The  partitions  of  T  are  illustrated  in  Figure  2.1. 

Assuming  the  kinematic  viscosity  to  be  constant  and  positive  (v  >  0),  the  following 
fluxes  are  defined: 


cr"  (u,  a,  p )  =  —  u  <g>  a  -  pi 

'advective'  flux 

(104) 

crd  (u)  =  2v\7su 

diffusive  flux 

(105) 

cr  (u,  a,  p)  =  crn  +  crd 

total  flux 

(106) 

cr“  =  n-  cra 

(107) 

crd  =  n-  crd 

(108) 

crn  =  n  ■  cr 

(109) 

where  p  is  the  pressure  (divided  by  the  density)  and  a  is  the  advective  velocity 

Solution  of  the  steady  Navier-Stokes  equations  involves:  find  u  =  u  (x) ,  p  =  p  (x) 

Q  such  that: 

Vx  G 

-V  •  cr{u,u,p)  =  f 

in  Q 

(110) 

V  •  u  =  0 

in  Q 

(111) 

u=  g 

°n  rg 

(112) 

—unu  —  pn  +  crd  (w)  =  h 

on  T;, 

(113) 

where  /  :  Q  — >  \Rd,  g  :  i\,  — *•  \Rd  and  h  :  T/,  — >  Rrf  are  prescribed.  The  boundary  condition 
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on  Th  can  be  interpreted  as: 


where 


h 


h  on  r;; 

< 

h+  on  r+ 


(114) 


crn  (u,  p)  =  h 

momentum  flux  boundary  condition 

(115) 

pn  +  cr^  (w,  p)  =  h+ 

traction  boundary  condition. 

(116) 

For  the  variational  formulation  of  the  Navier-Stokes  equations,  the  following  function 
spaces  are  required: 


j«e(wl(n))  \u=8  °n  T?  | 

(117) 

e  ^H1  (Q)^  \  w  =  0  on  r^| 

(118) 

{P  £  L2  (O)}. 

(119) 

The  variational  problem  consists  of:  find  u  G  <S,  p  G  V  such  that: 

B  ({w,q};{u,p},u)  =  L(w)  Vw  G  V,  G  V  (120) 


where 


B  ({w,q};{n,p},a)  =  (Vw,  a  (u,  a,  p))n  +  (w,a+u)Th  -  {Vq,u)Ci  +  (q,u„)T  (121) 


and 


L{w)  =  (w,f)n  +  (w,h)Th. 


(122) 
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Remark 


If  q  and  u  are  sufficiently  smooth,  -  (Vq,  u)n  +  (q,  un)T  =  (q,  V  •  «)n.  It  is  assumed  that 
this  is  the  case.  For  the  finite  element  case,  differentiability  and  continuity  will  be  ad¬ 
dressed  more  precisely  □ 

Consistency  of  the  weak  formulation  with  the  strong  form,  namely  equations  (110)  to 
(113),  can  be  proven  by  integrating  equation  (120)  by  parts,  leading  to: 

0  =  B  ({ w,q}-,{u,p},u )  —  L  (w) 

(123) 

=  -  (w,  V  •  cr  (u,u,  p)  +  /)n  +  ( w,  — un  u  —  pn  +  a „  ( u )  -  h )  +  (q,  V  •  u)n  . 

V  J  T/, 

Stability  requires  that  the  form  in  equation  (121)  be  coercive  in  the  following  sense.  Setting 
{u,  p}  =  {w,  q}  and  a  =  u  in  equation  (121), 

B  ({ w,q };  {w,q},u)  =  (Vw,  —w  ®  u  -  ql  +  2vVsw)Cii  +  (w,u%w)T  -  (V^,  w)n 

1  1 

=  -  (w,  (V  •  u )  w)n  —  -  ( w/unw)Ih  +2v||  Vsz^||q  +  {w,u^w)Yh 

1  1 

=  2  ( w >  (V  '  u)  w)n  +  2  lM”|w)r(i  +2v||Vsm||^. 

(124) 


This  can  be  restated  as: 

1  1 

B  ({ w,q};{w,q},u )  =  -  (w,  (V  •  u)  w)n  +  -|||w)!|1/2w||rji  +2i/||Vsm||^  (125) 

which,  for  the  case  of  a  divergence-free  flow  field  (V  •  u  =  0),  guarantees  stability  of 
the  field  w.  However,  analogous  to  the  advection-diffusion  equation,  the  problem  is  not 
necessarily  stable  if  the  flow  field  is  not  divergence-free. 
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3.2  Conservation  of  mass  and  momentum 


Conservation  of  momentum  is  examined  for  the  general  case  T,s>  Analogous  to  the 

advection-diffusion  case,  this  requires  the  definition  of  a  space  W  such  that 

(h1(C1)Y  =  V©W.  (126) 

The  problem  then  involves:  find  u  G  S,  p  €  V  and  hg  €  TV,  where  hg  is  the  momentum 
flux  across  Tg/  such  that 

B  ({zv,q};{u,p},u)  =  L  (zv)  +  (zv,hg)T  VmG^H1(Q)^  ,\/qeV.  (127) 

Setting  {w,q}  =  {e,,0\,  where  e,  is  the  canonical  unit  basis  vector,  in  the  above  equation. 


0  =  B  ({eit 0},  {«,  p},«)  -  (c,-,/)n  -  («//*») r„  “  (cA)r? 

=  /  u^ej-udT—  /  ei  fdCl—  /  ej  hdT—  /  e,  •  hg  dY. 
JTu  J  Cl  JT  U  JTq 


(128) 


This  can  be  rewritten  as 


0  =  /  ej-fdCl+  f  e,  h  dT  +  j  e,  ■  (— +  /i+)  dT  +  f  ej-hgdT  (129) 
kJ  o  v  r  ^  */  r  ^  r  ^ 

which  proves  that  momentum  is  conserved.  Setting  =  {0, 1}, 


0  =  B  ({0,1},  {«,  p},  w)  -  (0,/)n  -  (0,fcg)r?  -  (0,/r)rji 


(130) 


which  shows  that  mass  is  also  conserved.  Note  that  in  proving  conservation,  no  reliance 
has  been  placed  on  the  flow  field  being  strictly  divergence-free. 
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3.3  Conservation  for  the  advective  form  of  the  Navier- Stokes  equations 

In  practice,  the  conservation  form  of  the  Navier-Stokes  equations  is  rarely  used  in  Galerkin 
computations.  Application  of  the  conservation  form  results  in  spurious  oscillations  due 
to  the  V  •  u  term  in  equation  (124),  which  acts  as  a  distribution  of  sinks  and  sources.  To 
achieve  accuracy,  the  advective  form  of  the  Navier-Stokes  equations  has  been  found  much 
superior  and  is  usually  adopted.  Using  the  identity: 

V  •  (u  <g>  a)  =  a  ■  Vm  +  u  (V  •  a)  (131) 

the  conservation  and  advective  forms  of  the  Navier-Stokes  equations  are  interchangeable 
for  an  advective  field  a  which  is  divergence-free  everywhere  in  Q. 

Considering  the  advective  term  in  the  Navier-Stokes  equations  in  isolation, 

—  (Vzv,u  <g>  a)fi  =  (zv,  a  ■  Vw)n  +  (w,  (V  •  a)  «)n  —  ( w,anu)T .  (132) 

Inserting  equation  (132)  into  equation  (121),  and  assuming  that  a  is  divergence-free, 

Ba  ({w,q};  {u,  p},a)  =  (w,  a  •  Vw)n  +  ~V I  +  ^/Vsm)0  —  ( w,a~u)Tf 

-  (■ w,anu)Tg  -  (Vq,u)n  +  (q,un) r  (133) 

which  is  the  advective  weak  form  of  the  Navier-Stokes  equations.  It  is  this  form  which  is 
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usually  utilised  in  finite  element  analysis.  Setting  {w,q}  =  [e,,  0 }, 


0  =  Ba  ({ei,0};{u,p},u)  -  L  (e,)  -  (et, hg)  r> 

=  (eu u  ■  Vu)n  -  {eirunu)rh  ~  {™>unu)Tg  -  (eirf)n  -  (eirh)Yh  -  {eifhg)Tg 
=  -  {eif  (V  •  u)  u)n  +  (eifunu)T  -  (eifu~u)Th  -  {eirunu)Tg  -  (eirf)n  -  (eirh) Yh  -  (, eirhg)T 

=  -  (««/  (v  '  «)  M)n  +  (««/««  «)r4  “  («i//)n  -  («i/fc)rfc  -  • 

(134) 

Equation  (134)  implies  that  for  momentum  to  be  conserved  when  using  the  advective 
form  of  the  Navier-Stokes  equations,  it  is  necessary  that: 

I  e,  •  (V  •  u)  u  dd  =  0  (135) 

J  n 

which  resembles  the  weak  form  of  the  continuity  equation,  which  requires  that: 


in 


q  ■  (V  •  u)  dCl  =  0  V  q  £  V. 


Hence,  satisfaction  of  weak  continuity  and  requiring  that: 


(136) 


V  D  H1  (O) 


(137) 


guarantees  that  momentum  will  be  conserved  in  the  advection  form  of  the  Navier-Stokes 
equations. 

The  difficulty  that  arises  when  adopting  a  Galerkin  approach  by  replacing  the  relevant 
infinite-dimensional  function  spaces  with  finite  dimensional  counterparts  is  that  satisfac- 
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tion  of  the  inf-sup  stability  condition  will,  in  general,  require  that  [12]: 

Vh  yh  (138) 

where  in  a  finite  element  context  yh  is  the  span  of  the  individual  velocity  shape  func¬ 
tion  components  associated  with  all  nodes.  Satisfaction  of  the  inf-sup  condition  will  in 
general  preclude  momentum  conservation.  This  is  the  case  for  Galerkin  methods.  For 
stabilised  Galerkin  methods,  equation  (138)  is  often  satisfied,  but  the  continuity  equation 
is  stabilised,  hence  a  form  other  than  (136)  is  employed.  This  in  fact  is  a  key  advantage,  as 
will  become  clear. 

3.4  Conservation  in  Galerkin  and  stabilised  Galerkin  formulations  of  the 
Navier-Stokes  equations  -  multiscale  approach 

Conservation  can  be  restored  to  stabilised  Galerkin  methods  in  a  relatively  simple  manner. 
Consider  again  a  multiscale  decomposition  of  the  advective  flow  field,  a  =  a  +  a',  in 
which  a  is  the  coarse  scale  component  and  a '  is  the  fine  scale  component.  As  before, 
a'  =  0  on  T.  Inserting  the  multiscale  decomposition  into  equation  (133), 

Fa,ms  {{w,q};{u,p},a)  =  Ba  ({ w,q};{u,p},a )  +  (w,a'  •  Vw)n  (139) 

A  Galerkin  procedure  is  now  followed  in  which  Sh  C  5,  Vh  C  V  and  Vh  C  V.  It  is 
assumed  that  a  is  uh  G  Sh ,  and  a'  is  u' ,  which  will  be  defined  shortly.  Similar  to  before,  a 
space  W,!  is  defined,  such  that 

=  Vh  ©  Wh  (140) 

where,  by  assumption,  Zh  =  (yh)d.  A  Galerkin  procedure,  including  a  multiscale  de¬ 
composition  of  the  advective  velocity,  involves:  find  uh  G  Sh,  ph  G  Vh  and  hg  G  Wh  such 
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that: 


B„  [[wh,qh}-,{uh,ph},uh^j  +  (wh,u'  ■  Vw^ 

=  L  (wh^j  +  (wh,h^j^  V  wh  €  Zh,  qh  €  Vh.  (141) 

Inserting  {wh,qh}  =  {e/,0}  into  equation  (141)  leads  to: 

0  =Ba  ^{e,,0};  {uh, ph},u^j  +  ( eiru ’  ■  Vm?!)^  -  L  (e/)  -  (e,,/?^ 

=  -  (eo  (V  •  w,!)  I/,!)n  +  (ei,u'  •  Vl/?‘)n  -  (e,-,/i_)r-  (142) 

-  (e/,  —u„uh  +  +  -  -  (c«,/)n- 

Hence,  momentum  is  conserved  if: 


((v-uh)  uh-u'-Vuh) 


dCl  =  0, 


and  this  requirement  is  satisfied  if: 


(143) 


J  (V'  (V  •  uh J  -  u!  ■  Vqhj  dd  =  0  V  qh  €  (144) 

subject  to  the  condition  that  Vh  D  yh.  Hence,  if  the  modified  continuity  equation  in  (144) 
is  satisfied  and  Vh  2>  yh,  the  formulation  in  equation  (141)  conserves  momentum. 

The  missing  component  to  this  point  is  the  definition  of  the  fine  scale  velocity.  This 
requires  the  supposition  of  a  model.  In  the  same  vein  as  for  the  advection-diffusion  prob¬ 
lem,  an  approximation  of  the  fine  scale  velocity  on  element  interiors  is  introduced  [2]: 


u'  =  — T pV 


(145) 
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where  r  is  the  residual  of  the  incompressible  Navier-Stokes  equation: 


r  =  uh  ■  Vuh  +  Vph  -  2vAuh  -  /.  (146) 

Inserting  the  expression  for  the  model  for  the  fine  scales  (145)  into  the  modified  continuity 
equation  (144)  yields: 

J  (V  (V  •  uh)  +  Vqh  ■  r -pr)  dCl  =  0  V  qh  £  Vh  (147) 

which  is  precisely  the  pressure-stabilised  continuity  equation.  Therefore,  a  momentum 
conserving,  pressure-stabilised  Galerkin  formulation  in  advection  format  involves:  find 
uh  G  Sh,  ph  G  'P,!  and  hg  G  >V/z  such  that: 

Bn  ([wh,qh};{uh,ph},uh^  +  (vqh, -  (V!,  rpr  •  Vh?'^ 

=  L  (w1')  +  (ivh,hg^  Vwh  £  Zh,Vqh  £Vh.  (148) 

While  the  above  formulation  is  conservative,  it  is  not  stabilised  for  advection  dominated 
flows.  The  following  stabilised  Galerkin  problem  is  both  stable  and  momentum  conserv¬ 
ing:  find  uh  £  Sh,  ph  £  Vh  and  hhg  £  Wh 

Bn  {{wh  ,qh}-,  {uh  ,ph},uh^  +  (vqh,rprSj^+  (cmwh,Tmr) 

+  ( wh,u '  ■  +  [L'mwh,x'  (u!  ■  Vnh))^ 

-  ( wh,h =  L  (V)  Vm/!  €  Z/!,  Vq,!  g  'P7'  (149) 

in  which  u'  =  — "lyr,  and 


(')  =  w'  •  V  (•)  =  — rpr  •  V  (•) . 


(150) 
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The  operator  Cm  defines  the  particular  stabilisation  method. 


uh  ■  Vwh 

(SUPG) 

uh  ■  Vwh  -  2vA wh 

(GLS) 

uh  ■  V w h  +  2vAivh 

(MS). 

(151) 


As  for  the  advection-diffusion  case,  the  term  introduced  to  provide  conservation  (the  first 
term  in  the  box)  is  an  advection  term  and  is  hence  stabilised  (second  term  in  the  box).  The 
second  term  in  the  box  can  be  interpreted  in  the  same  fashion  as  equation  (67).  Appropri¬ 
ate  expressions  for  rp,  rm  and  x'  can  be  found  in  Taylor  et  al.  [13],  in  which  it  is  assumed 
that  xm  =  x p. 


3.5  Time  dependent  case 

Considering  the  space-time  domains  defined  in  equations  (68)-(74),  solution  of  the  un¬ 
steady  Navier-Stokes  equations  involves:  find  u  =  u  (x,  t) ,  p  =  p(x,t )  Vx  G  O,  Vf  G  [0,  T] 
such  that: 


V  •  cr  (■ u ,  u,  p)  =  f 

in  Q 

(152) 

u  (x,  0)  =  «o 

in  O 

(153) 

V  •  u  =  0 

in  Q 

(154) 

u=  g 

on  Pg 

(155) 

—  pn  +  ai  (w)  =  h 

on  Ph 

(156) 

where  /  :  Q  — >  Rd,  g  :  Pg  — *  Rd  and  h  :  Ph  — >  Rrf  are  prescribed. 

Considering  the  space-time  'slabs'  defined  in  equations  (79)-(85),  the  pertinent  function 
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spaces  are: 


{“€( 

H1  (Q”)y?  w  =  g  on  Pg | 

(157) 

|m  £  | 

^H1(Q',))'f  =  0  on  Pg | 

(158) 

{p£L 

2(Q")} 

(159) 

\d 

m) 

=  V'!  ©  W" 

(160) 

and  relevant  subspaces  are: 


(Sn)h  C  Sn  (161) 

( Vn)h  C  V”  (162) 

( Vn)h  C  Tn  (163) 

(« Sn)h  C  ( Zn)h  =  ( Vn)h  ©  (Wn)h  .  (164) 


Solution  of  the  unsteady  problem  requires:  find  uh  £  (Sn)h ,  ph  £  (Vn)n  such  that 


Q" 


+  (  uAu''  •  Vnh 


Q" 


+  (  Vie/!,  —pi  +  2vVuh  j  ^  —  ^m?!/  wjj  u  j  ^  u\uh ^  ^ 


:  -  Ky,)0. 


+  (^w«)p„  -  (wh'hg)pn  +  •  V«*)q.  +  (^wfc,T'  («'  •  Vm/!))q„ 


( Vqh,u n  +  Cmwhfrmrj  ~n  =  L  (wh^j 


Vw11  £  (Zn)h,  Vqh  £  (7>n)\  n  =  0, 1, . . . ,  N  —  1  (165) 


where 


(166) 
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and  the  fine  scale  velocity  is  given  by  u'  =  —  Tpr  and  C!m  is  defined  by  equation  (150).  The 
operator  Cm  for  the  unsteady  case  is  defined  by: 


Cmwh 


dwh 

dt 

dwh 

dt 

dwh 
„  dt 


+  uh  ■  Vwh 
+  uh  -Vwh  —  vA wh 
+  uh  -Vivh+  vA wh 


(SUPG) 

(GLS) 

(MS). 


(167) 


The  boxed  terms  are  the  same  as  for  the  steady-state  case  (see  equation  (149)).  Considering 
first  {wh,qh}  =  {0, qh }  and  inserting  the  model  u'  =  —  rpr, 


0  =  -  (v^V)  +  (<;''.  "S)p„  +  K'.T„r)Sii 

(160) 

which  is  equation  (147),  the  weak  form  of  the  stabilised  continuity  equation. 

Remark 

Conservation  of  mass  may  be  obtained  from  the  first  line  of  equation  (168)  by  selecting 
qh  =  1.  The  integration-by-parts  utilised  in  the  second  line  of  equation  (168)  assumes 
that  qh  is  continuous  across  element  interfaces  and  exact  quadrature  is  performed  (uh  is 
always  assumed  to  be  continuous).  If  qh  is  discontinuous  across  element  interfaces,  the 
first  line  of  equation  (168)  requires  additional  pressure  jump  terms  on  element  interfaces, 
whereas  the  second  line  is  still  valid.  This  follows  by  performing  integration-by-parts 
element-wise: 


q\  V  •  uh 


Qn 


( Vqh,i 


(169) 


+  (  (??!/  un  )  + 


qhn 
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where 


qhl_n+  +  c\h_n __ 

(qh+-qh-)  n+ 

-  (?+-?-)  n- 


in  which  [[•]]  is  the  jump  operator  across  element  interfaces  f. 


(170) 


r  =  (Jre\r,  p”  =  rex]tn/tn+1[,  p"  =  fx]  tn,tn+1[  (i7i) 

e 

and  the  ±  designations  are  arbitrary  labels  of  quantities  associated  with  two  different 
elements  sharing  an  interface.  Note  that  equation  (170)  is  invariant  under  reversal  of  the 
±  labels.  When  inexact  element  quadrature  is  used  for  the  continuous  qh  case,  the  form  in 
the  first  line  of  equation  (168)  is  preferable  because  conservation  is  unaffected.  This  is  not 
necessarily  the  case  with  the  form  on  the  second  line  of  (168).  Perhaps  because  most  code 
developers  are  unaware  of  this  fact,  the  second-line  form  is  often  employed  in  practice.  □ 
Following  the  now  familiar  procedure  of  setting  {ivh,qh}  =  {<?/,  0 } ,  and  inserting  this 
into  the  unsteady  Navier-Stokes  equation. 
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Using  the  result  in  equation  (168),  the  above  equation  is  equivalent  to: 


+  [  erh~dP+  I  ei-hhrdP+  [  et  ■  f  dQ  (173) 

Jp-n  JP'l  c  Jq " 

which  proves  that  for  the  time-dependent  case  conservation  of  momentum  is  attained  on 
space-time  slabs. 

4  Conclusions 

The  conservation  properties  of  Galerkin  and  some  common  stabilised  Galerkin  schemes 
for  advection-diffusion  and  Navier-Stokes  equations  have  been  examined.  It  was  shown 
why  advection-diffusion  and  Navier-Stokes  problems  based  on  advective  weak  forms, 
computed  from  Galerkin  and  stabilised  Galerkin  finite  element  methods,  are  not  gener¬ 
ally  conservative.  It  was  however  shown  that,  through  multiscale  considerations,  this  de¬ 
ficiency  can  be  easily  overcome  for  stabilised  methods  by  the  addition  of  a  residual  term 
to  the  weak  equation  which  ensures  global  conservation.  Surprisingly,  conservation  can 
be  restored  more  easily  to  stabilised  Galerkin  methods  than  traditional  Galerkin  methods. 
The  reason  for  this  is  that  conservation  is  in  conflict  with  the  inf-sup  condition  associated 
with  incompressibility  in  the  Galerkin  formulation,  whereas  in  stabilised  formulations  the 
inf-sup  condition  is  circumvented. 

Stabilised,  conservative  formulations  have  been  described  and  presented  for:  steady 
advection-diffusion  equation  (equation  (64));  unsteady  advection-diffusion  equation  (equa¬ 
tion  (92));  steady  incompressible  Navier-Stokes  equations  (equation  (149));  and  unsteady 
incompressible  Navier-Stokes  equations  (equation  (165)).  The  techniques  described  herein 
have  been  extensively  tested,  and  verified  in  several  academic  and  commercial  flow  codes. 
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including  Spectrum  [13]  and  Acusolve  [14]. 

Similar  issues  are  being  investigated  in  the  discontinuous  Galerkin  method  literature. 
The  interested  reader  is  referred  to  Cockburn  et  al.  [15, 16]. 
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